arXiv:1509.05377vl [cs.CG] 17 Sep 2015 


Computing the Rectilinear Center of Uncertain 
Points in the Plane 


Haitao Wang and Jingru Zhang 

Department of Computer Science 
Utah State University, Logan, UT 84322, USA 
haitao.wangSusu.edu,jingruzhangSaggiemail.usu.edu 


Abstract. In this paper, we consider the rectilinear one-center problem 
on uncertain points in the plane. In this problem, we are given a set V 
of n (weighted) uncertain points in the plane and each uncertain point 
has m possible locations each associated with a probability for the point 
appearing at that location. The goal is to find a point q* in the plane 
which minimizes the maximum expected rectilinear distance from q* to 
all uncertain points of V, and q* is called a rectilinear center. We present 
an algorithm that solves the problem in 0{mn) time. Since the input 
size of the problem is &(mn), our algorithm is optimal. 


1 Introduction 

In the real world, data is inherently uncertain due to many facts, such as 
the measurement inaccuracy, sampling discrepancy, resource limitation, and so 
on. A large amount of work has recently been done on uncertain data, e.g., 
|ll2l;il9ll2ll:ilt7m^ . In this paper, we study the one-center problem on uncer¬ 
tain points in the plane with respect to the rectilinear distance. 

Let V = {Pi, P 2 , ■ ■ ■ ,Pn} be a set of n uncertain points in the plane, where 
each uncertain point Pi € V has m possible locations pn , Pi 2 ,■■■ , Pim and for 
each 1 < j < m, pij is associated with a probability > 0 for Pi being at pij 
(which is independent of other locations). 

For any (deterministic) point p in the plane, we use Xp and yp to denote the x- 
and ^-coordinates of p, respectively. For any two points p and q, we use d(p, q) to 
denote the rectilinear distance between p and q, i.e., d{p,q) = \xp — Xq\ + \yp — yq\. 

Consider a point q in the plane. For any uncertain point Pi S P, the expected 
rectilinear distance between q and Pi is defined as 

m 

Ed{P^, g) = X fb ■ d{Ptj,q)- 
i=i 

Let Edmax(g) = rnaxp^g-p Ed(Pi, g). A point q* is called a rectilinear center 
of P if it minimizes the value Edniax(g*) among all points in the plane. Our goal 
is to compute g*. Note that such a point g* may not be unique, in which case 
we let g* denote an arbitrary such point. 




We assume that for each uncertain point Pi oi V, its m locations are given 
in two sorted lists, one by ^-coordinates and the other by y-coordinates. To the 
best of our knowledge, this problem has not been studied before. In this paper, 
we present an 0{mn) time algorithm. Since the input size of the problem is 
0{nm), our algorithm essentially runs in linear time, which is optimal. 

Further, our algorithm is applicable to the weighted version of this problem 
in which each Pi € V has a weight Wi > 0 and the weighted expected distance^ 
i.e., Wi ■ Ed(Pi, q), is considered. To solve the weighted version, we can first reduce 
it to the unweighted version by changing each to Wi ■ fij for all 1 < f < n 
and 1 < j < m, and then apply our algorithm for the unweighted version. The 
running time is still 0{mn). 

1.1 Related Work 

The problem of finding one-center among uncertain points on a line has been 
considered in our previous work m, where an 0{nm) time algorithm was given. 
An algorithm for computing k centers for general k was also given in |21) with 
the running time 0(mn log mn-|-n log n log fc). In fact, in [21] we considered the 
fc-center problem under a more general uncertain model where each uncertain 
point can appear in m intervals. We also studied the one-center problem for 
uncertain points on tree networks in |20j . where a linear-time algorithm was 
proposed. 

There is also a lot of other work on facility location problems for uncertain 
data. For instances, Cormode and McGregor [7] proved that the /c-center problem 
on uncertain points each associated with multiple locations in high-dimension 
space is NP-hard and gave approximation algorithms for different problem mod¬ 
els. Foul m considered the Euclidean one-center problem on uncertain points 
each of which has a uniform distribution in a given rectangle in the plane, de 
Berg, et al. [8] studied the Euclidean 2-center problem for a set of moving points 
in the plane (the moving points can be considered uncertain). 

The /c-center problems on deterministic points are classical problems and have 
been studied extensively. When all points are in the plane, the problems on most 
distance metrics are NP-hard m- However, some special cases can be solved in 
polynomial time, e.g., the one-center problem [H], the two-center problem [B], 
the rectilinear three-center problem m, the line-constrained /c-center problems 
(where all centers are restricted to be on a given line in the plane) |5II4II9) . 

1.2 Our Techniques 

Consider any uncertain point Pi G P and any (deterministic) point q in the 
plane We first show that Ed(Pi, q) is a convex piecewise linear function with 
respect to g G R^. More specifically, if we extend a horizontal line and a vertical 
line from each location of Pi, these lines partition the plane into a grid Gi of 
(to -I- 1) X (m -I- 1) cells. Then, Ed(Pi, g) is a linear function (in both the x- and 
^-coordinates of q) in each cell of Gi. In other words, Ed{Pi,q) defines a plane 
surface patch in 3D on each cell of Gi. Then, finding q* G R^ is equivalent to 
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finding a lowest point p* in the upper envelope of the n graphs in 3D defined by 
Ed(Pi, q) for all Pi GV (specifically, q* is the projection of p* onto the xy-plane). 

The problem of finding p*, which may be interesting in its own right, can 
be solved in 0(nm?') time by the linear-time algorithm for the 3D linear pro¬ 
gramming (LP) problem |15) . Indeed, for a plane surface patch, we call the plane 
containing it the supporting plane. Let T-L be the set of the supporting planes of 
the surface patches of the functions Ed(Pi, q) for all Pi G V. Since each function 
Ed(Pi, q) is convex, p* is also a lowest point in the upper envelope of the planes 
of T-L. Thus, finding p* is a LP problem in and can be solved in 0(|H|) time 
[T5] . Note that \'H\ — 0{nw?) since each grid Gi has (m + 1)^ cells. 

We give an 0{mn) time algorithm without computing the functions Ed(Pi, q) 
explicitly. We use a prune-and-search technique that can be considered as an 
extension of Megiddo’s technique for the 3D LP problem [15]. In each recursive 
step, we prune at least n/32 uncertain points from P in linear time. In this way, 
q* can be found after O(logn) recursive steps. 

Unlike Megiddo’s algorithm each recursive step of our algorithm itself 
is a recursive algorithm of 0(\ogm) recursive steps. Therefore, our algorithm 
has O(logn) “outer” recursive steps and each outer recursive step has O(logm) 
“inner” recursive steps. In each outer recursive step, we maintain a rectangle R 
that always contains q* in the xp-plane. Initially, R is the entire plane. Each 
inner recursive step shrinks R with the help of a decision algorithm. The key 
idea is that after 0(log m) steps, R is so small that there is a set V* of at least 
n/2 uncertain points such that R is contained inside a single cell of the grid Gi 
of each uncertain point Pi of V* (i.e., R does not intersect the extension lines 
from the locations of Pi). At this point, with the help of our decision algorithm, 
we can use a pruning procedure similar to Megiddo’s algorithm m to prune 
at least |P*|/16 > n/32 uncertain points of V*. Each outer recursive step is 
carefully implemented so that it takes only linear time. 

In particular, our decision algorithm is for the following decision problem. Let 
i? be a rectangle in the plane and R contains q* (but the exact location of q* is 
unknown). Given an arbitrary line I that intersects i?, the decision problem is to 
determine which side of I contains q*. Megiddo’s technique m gave an algorithm 
that can solve our decision problem in 0(m?n) time. We give a decision algorithm 
of 0{mn) time. In fact, in order to achieve the overall 0{mn) time for computing 
g*, our decision algorithm has the following performance. For each 1 < * < n, 
let Oi and bi be the number of columns and rows of the grid Gi intersecting R, 
respectively. Our decision algorithm runs in 0(^/^j^(ai -I- bi)) time. 

The rest of the paper is organized as follows. In Section [5] we introduce some 
observations. In Section[3l we present our decision algorithm. Section |4| gives the 
overall algorithm for computing the rectilinear center q*. Section [S] concludes. 

2 Observations 

Let p be a point in the plane M^. The vertical line and the horizontal line through 
p partition the plane into four (unbounded) rectangles. Consider another point 
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9 G We consider d{p,q) as a function of g G For each of the above 
rectangle R, d(j), g) on g G i? is a linear function in both the x- and ^-coordinates 
of q, and thus d{p, q) onq € R defines a plane surface patch in R^. Further, d{p, q) 
on g G R^ is a convex piecewise linear function. 

For ease of exposition, we make a general position assumption that no two 
locations of the uncertain points of V have the same x- or y-coordinate. 

Consider an uncertain point Pi oi V. We extend a horizontal line and a 
vertical line through each location of Pi to obtain a grid, denoted by Gi, which 
has (m + 1) X (m -|- 1) cells (and each cell is a rectangle). According to the 
above discussion, for each location pij of P, the function d{pij,q) of g in each 
cell of Gi is linear and defines a plane surface patch in R^. Therefore, if we 
consider Ed(Pi,g) as a function of g, since Ed{Pi,q) is the sum of fij • d(pij,q) 
for all 1 < j < TO, Ed(Pi,g) of g in each cell of Gi is also linear and defines a 
plane surface patch in R^. Further, since each d{pij,q) for g G R^ is convex, the 
function Ed(Pi,g), as the sum of convex functions, is also convex. 

In the following, since Ed{Pi,q) is normally considered as function of g, for 
convenience, we will use Edi(a;,?/) to denote it for g = {x,y) G R^. 

The above discussion leads to the following observation. 

Observation 1 For each uncertain point Pi G V, the function Edi(x, y) is con¬ 
vex piecewise linear. More specifically, Edi(x,y) on each cell of the grid Gi is 
linear and defines a plane surface patch in R^ (e-Q-, see Fig. Qp. 

Consider the function Edi{x,y) of any Pi G V. Clearly, the complexity of 
Edi{x,y) is 0{wf). However, since Edi{x,y) on each cell G of Gi is a plane 
surface patch in R^, Edi(x, y) on G is of constant complexity. We use Edi{x, y, G) 
to denote the linear function of Edi(x, y) on C. Note that Edi(x, y, C) is also the 
function of the supporting plane of the surface patch of Edi(x, y) on C. 

As discussed in Section o our algorithm will not compute the function 
Edi{x,y) explicitly. Instead, we will compute it implicitly. More specifically, 
we will do some preprocessing such that given any cell C oi Gi, the function 
Edi{x,y,G) can be determined efficiently. We first introduce some notation. 

Let Xi = {xii,Xi 2 , ■ ■ ■ ,Xim} be the set of the x-coordinates of all locations 
of Pi sorted in ascending order. Let Yi = {yn, ya, - • • , Vim} be the set of their 
^-coordinates in ascending order. Note that Xi and Yi can be obtained in 0{m) 
time from the input (recall that the locations of Pi are given in two sorted lists 
in the input). For convenience of discussion, we let Xio = —oo, and let Xi also 
include Xio. Similarly, let ym = —oo, and let Yi also include yio. Note that due 
to our general position assumption, the values in Xi (resp., Yi) are distinct. 

For any value z, we refer to the largest value in Xi that is smaller or equal 
to z the predecessor of z in Xi, and we use Iz{Xi) to denote the index of the 
predecessor. Similarly, IzfYi) is the index of the predecessor of z in Yi. 

Consider any point g in the plane. The predecessor of the x-coordinate of g 
in Xi is also called the predecessor of q in Xi. Similarly, the predecessor of the 
j/-coordinate of g in Yi is also called the predecessor of q in Yi. We use Iq(Xi) 
and Iq{Yi) to denote their indices, respectively. 
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Figure 1. Illustrating the function EAi(x,y) of an uncertain point Pi with m = 4. 

Consider any cell C of the grid Gi. For convenience of discussion, we assume 
C contains its left and bottom sides, but does not contain its top and right sides. 
In this way, any point in the plane is contained in one and only one cell of Gi. 
Further, all points of G have the same predecessor in Xi and also have the same 
predecessor in Yi. This allows us to define the predecessor of G in Xi as the 
predecessor of any point in Xi, and we use Ic{Xi) to denote the index of the 
predecessor. We define IciXi) similarly. We have the following lemma. 

Lemma 1. For any uncertain point Pi € V, after 0{m) time preprocessing, for 
any cell G of the grid Gi, if Ic{Xi) and Ic(Yi) o,re known, then the function 
Ed i{x,y,G) can be computed in constant time. 

Proof. For each location p € Pi, let Xp and pq be the x- and y-coordinates of p, 
respectively, and let fp be the probability associated with p. 

For any point q = {x,y) in K^, recall that the expected distance function 
Ed^{x,y) = T,pePi fp-d{p,q) = T,pePi fp-i\^p-^\ + \yp-y\)- Therefore, we can 
write Edi{x,y) = J2pePi fp'\^p-^\+J2pePi fp-\yp-y\ - In Hi® following, we first 
discuss how to compute J^p^Pi fp ' “ ^1 nnd the case for J^p^Pi fp ' \yp ~ y\ 

is very similar. 

Let S'! denote the set of all locations of Pi whose x-coordinates are smaller 
than or equal to x, i.e., the ^-coordinate of q. Let S 2 = Pi \ Si. Then, we have 
the following: 



psPi 


P&Si 


P&S 2 



( 1 ) 


pGSi p£52 p^S2 


P^S2 


— X ■ (2 ■ fp /p) 2 fp ■ Xp + fp ■ Xp. 


P^Si p^Pi p^Si p^Pi 


p^pi 
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Thus, in order to compute fp' l^p it is sufficient to know the four 

values EpGSi fp’ SpGPi fp’ SpeSi fp ' SpePi fp ' ^p- this end, we do 

the following preprocessing. 

First, we compute X^peP^ fp J2p&Pi fp ' ^p) which can be done in 0(m) 
time. Second, recall that Xi = {xio,Xii,... ,Xim} maintains the ^-coordinates 
of the locations of Pi sorted in ascending order. Note that given any index j 
with 1 < j < TO, we can access the information of the location of Pi whose 
x-coordinate is Xij in constant time, and this can be done by linking each Xij 
to the corresponding location of Pi when we create the list Xi from the input. 
For each j with 1 < j < n, we let f{j) be the probability associated with the 
location of Pi whose a:-coordinate is Xij. 

In the preprocessing, we compute two arrays A[0 • • • to] and i?[0 • • • to]. Specif¬ 
ically, for each 1 < j < n, A[j] = /(k) and B[j] = fik) ■ x,k- For 

j = 0, we let A[0] = i3[0] = 0. As discussed above, since we can access /(j) in 
constant time for any 1 < j < to, the two arrays A and B can be computed in 
0{m) time. 

Let t = Iq{Xi), i.e., the index of the predecessor of q in Xi. Note that 
t € [0, to]. To compute J2^Pi /p' I^p “^1 > easy observation is that J^p^Si fp 
exactly equal to A[t\ and fp'Xp is exactly equal to B[t\. Therefore, with the 

above preprocessing, if t is known, according to Equation ([T]), J2p&Pi fp ' l^p ~^\ 
can be computed in 0(1) time. 

The above shows that with 0{m) time preprocessing, given Iq{Xi), we can 
compute the function X^pGP fp ’ I^p “ of a: at g = (x, y) in constant time. 

In a similar way, with Oim) time preprocessing, given IqiYi), we can compute 
the function X^peP^ fp ' IVp ~ J/l of y at p = (x, y) in constant time. 

Let q be any point in the cell C. Hence, Iq{Xi) = Ic{Xi) and IqiYi) = 
IciYi). Further, the function Edi{x,y) on q = (x,y) G O is exactly the function 
Edi(x, y, C). Therefore, with 0 (to) time preprocessing, given Ic{Xi) and IciYi), 
we can compute the function Edi(x,y, O) in constant time. 

The lemma thus follows. □ 

Due to Lemma [TJ we have the following corollary. 

Corollary 1. For each uncertain point Pi G V, after Oim) time preprocessing, 
given any point q in the plane, the expected distance EdiPi,q) can be computed 
in O(logTO) time. 

Proof. Given any point q G we can compute IqiXi) in O(log to) time by 
doing binary search on Xi. Similarly, we can compute IqiYi) in O(logTO) time. 
Let C be the cell containing q. Recall that IciXi) = IqiXi) and IciXi) = IqiXi)- 
Hence, by Lemma [1] we can compute the function Edi(x, y, C) in constant time. 
Then, Ed(Pi,q) is equal to EdiX, qy, C), where qx and qy are the x- and y- 
coordinates of q, respectively. Thus, after Edi(x, y, C) is known, Ed(Pi, q) can be 
computed in constant time. The corollary thus follows. □ 

Recall that Edinax(q) = maxp^g-p Ed(Pi,q) for any point q in the plane. For 
convenience, we use Edmax(a^j y) to represent Edniax(q) as a function of q = 
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(x, y) € Note that Edmax(a;, y) is the upper envelope of the functions Edi(x, y) 
for alH = 1, 2,..., n. Since each Edi{x,y) is convex on Edniax(a^, y) is also 
convex on Further, the rectilinear center q* corresponds to a lowest point p* 
on Edniax(a^j y)- Specifically, q* is the projection of p* on the xj/-plane. Therefore, 
computing q* is equivalent to computing a lowest point in the upper envelope 
of all functions Edi(a:, y) for all z = 1, 2,..., n. 

For each 1 < z < n, let Hi denote the set of supporting planes of all surface 
patches of the function Edi(a::, 2 /). Let H. = Since Edi{x,y) is convex, 

Edi(a;, y) is essentially the upper envelope of the planes in Hi. Hence, Edmaxix, y) 
is also the upper envelope of all planes in H. Therefore, as discussed in Section 
11.21 finding p* is essentially a 3D LP problem on H, which can be solved in 
0(1^1) time by Megiddo’s technique [15]. Since the size of each Hi is 0(rn?), 
|H| = 0{nm?). Therefore, applying the algorithm in [15] directly can solve the 
problem in 0{nrn?) time. In the following, we give an 0{nm) time algorithm. 

In the following paper, we assume we have done the preprocessing of Lemma 
[Tjfor each Pi G P, which takes 0{mn) time in total. 

3 The Decision Algorithm 

In this section, we present a decision algorithm that solves a decision problem, 
which is needed later in Section 21 We first introduce the decision problem. 

Let R = [xi,X 2 ;yi, 2 / 2 ] be an axis-parallel rectangle in the plane, where xi 
and X 2 are the cc-coordinates of the left and right sides of R, respectively, and 
yi and ?/2 are the ^-coordinates of the bottom and top sides of R, respectively. 
Suppose it is known that q* is in R (but the exact location of q* is not known). 
Let L be an arbitrary line that intersects the interior of R. The decision problem 
asks whether q* is on L, and if not, which side of L contains q*. We assume the 
two predecessor indices I^^^Xi) and Iyj^(Yi) are already known. 

For each I < i < n, let Oi = Ix.^{Xi)—Ix^ {Xi) + 1 and bi = Iy.^(Yi)—Iyj^{Yi) + l. 
In fact, Qi and bi are the numbers of columns and rows of Gi that intersect i?, 
respectively. Below, we give a decision algorithm that solves the decision problem 
in d-bi)) time. Note that 2n < + h) < 2{m -\- l)n. 

We first show that the decision problem can be solved in 0(X]r=i ' ^0 time 

by using the decision algorithm for the 3D LP problem m- Later we will reduce 
the running time to 0(X]r=i(®* + ^*)) f™®- 

Recall that p* is a lowest point in the upper envelope of the functions Edi(a;, y) 
for * = I, 2,..., n. Since q* is in R and each function Edi(a;, y) is convex, an easy 
observation is that p* is also a lowest point in the upper envelope of Edi{x,y) 
for i = 1,2,... ,n restricted on (x, y) G R. This implies that we only need to 
consider each function Edi(x,y) restricted on R. 

For each I < z < n, let Gi{R) be the set of cells of Gi that intersect R, 
and let Hi{R) be the set of supporting planes of the surface patches of Ed(Pi, q) 
defined on the cells of Gi{R). Let 'H{R) = By our above analysis, 

p* is a lowest point of the upper envelope of all planes in T-L{R). Note that 
\Hi{R)\ = m ■ bi for each I < z < n. Thus, \'H{R)\ = Then, we 
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can apply the decision algorithm in [TS] (Section 5.2) on to determine 

which side of L contains q* in 0{\'H{R)\) time. In order to explain our improved 
algorithm later, we sketch this algorithm below. 

We consider each plane of 'H{R) as a function of the points q on the cry-plane 
K.^. In the first step, the algorithm finds a point q' on L that minimizes the 
maximum value of all functions in R{R) restricted on the line q € L. This is 
essentially a 2D LP problem because each function of H(i?) restricted on L is 
a line, and thus the problem can be solved in 0{\'H{R)\) time [15]. Let (pqi be 
the set of functions of 'H{R) whose values at q' are equal to the above maximum 
value. The set (Pq' can be found in 0{\R{R)\) time after q' is computed. This 
finishes the first step, which takes 0{\R{R)\) = 0(X]r=i 

The second step solves another two instances of the 2D LP problem on 
the planes of <l>q', which takes 0{\<l>qi\) time. An easy upper bound for 
is X]r=i ^ close analysis can show that \<l>qi\ = 0{n). Indeed, for each 

1 < i < n, since the function Edi(cr,y) is convex, among all ai ■ bi planes in 
Hi{R), at most four of them are in Therefore, \'Pq' \ = 0{n). Hence, the sec¬ 
ond step runs in 0{n) time. Since in our problem there always exists a solution, 
according to m, the second step will either conclude that q' is q* or tell which 
side of L contains q*, which solves the decision problem. The algorithm takes 
Qi ■ bi) time in total, which is dominated by the first step. 

In the sequel, we reduce the running time of the above algorithm, in par¬ 
ticular, the first step, to 0{J27=ii^i + ^i))- Our goal is to compute q' and <l>q'. 
By the definition, q' is a lowest point in the upper envelope of all functions 
of R{R) restricted on the line L. Consider any uncertain point Pi € V. Let 
Hi{R,L) be the set of supporting planes of the surface patches defined on the 
cells of Gi{R) intersecting L. Observe that since Edi{x,y) is convex, the upper 
envelope of all the functions of Hi{R) restricted on L is exactly the upper enve¬ 
lope of the functions of Hi{R,L) restricted on L. Therefore, q' is also a lowest 
point in the upper envelope of the functions of L) restricted on L, where 
'H{R,L) = L). In other words, among all planes in R{R), only the 

planes of 'H{R,L) are relevant for determining q' . Thus, suppose 'H{R,L) has 
been computed; then q' can be computed based on the planes of R{R,L) in 
0{\'}i{R,L)\) time by the 2D LP algorithm [15]. After q' is computed, the set 
can also be determined in 0{\H{R,L)\) time. 

Note that L)\ = 0(X]r=i(®* + since for each 1 < f < n, \Hi{R, L)|, 

which is equal to the number of cells of Gi{R) intersecting L, is 0{ai + bi). 

It remains to compute 'H{R,L), i.e., compute Hi{R,L) for each 1 < z < n. 
Recall that R = [xi,X 2 \yi,y 2 ] and the two predecessor indices R^iXi) and 
Iyi(Yi) for each 1 < z < n are already known. The following lemma gives an 
0{ai + bi) algorithm to compute Hi{R,L). 

Lemma 2. For each 1 < z < zz, Hi{R,L) can be computed in 0{ai bi) time. 

Proof. Let Gi{R,L) be the set of cells of Gi{R) intersecting L. To compute 
the planes in Hi{R,L), it is sufficient to determine the plane surface patches 
of Edi{x,y) defined on the cells of Gi{R,L). By Lemma (U this amounts to 


determine the indices of the predecessors of these cells in Xi and Yi, respectively. 
In the following, we give an algorithm to compute the cells of Gi{R,L) and 
determine their predecessor indices in Xi and Yi, respectively, and the algorithm 
runs in 0(oi + bi) time. 

The main idea is that we first pick a particular point p on L (1 R and locate 
the cell of Gi{R) containing p (clearly this cell belongs to Gi{R,L)), and then 
starting from p, we traverse on L and Gi (i?) simultaneously to trace other cells 
of Gi{R, L) until we move out of R. The details are given below. 

We focus on the case where L has a positive slope. The other cases can be 
handled similarly. Recall that L intersects the interior of R. Let p be the leftmost 
intersection of L with the boundary of R. Hence, p is either on the left side or 
the bottom side of R. 

Let G be the cell of Gi that contains p. We first determine the two indices 
/p(W) and /p(r,). Note that IciX,) = Ip{X,) and Ic(Y^ = Ip{Y,). 

Since p € R, the index Ip(Xi) can be found in 0(ai) time by scanning the 
list Xi from the index Rj^iXi). Similarly, Ip{Yi) can be found in 0{bi) time by 
scanning the list Yi from the index ly^iYi). After Ic{Xi) = Ip{Xi) and IciXi) = 
Ip{Yi) are computed, by Lemma[I] the function Edi(a;, y, G) can be computed in 
constant time, and we add the function to Hi{R,L). 

Next, we move p on L rightwards. We will show that when p crosses the 
boundary of G, we can determine the new cell containing p and update the two 
indices Ip(Xi) and Ip{Yi) in constant time. This process continues until p moves 
out of R. Specifically, when p moves on L rightwards, p will cross the boundary 
of G either from the top side or the right side. 

First, we determine whether p will move out of R before p crosses the bound¬ 
ary of G. If yes, then we terminate the algorithm. Otherwise, we determine 
whether p moves out of G from its right side or left side. All above can be easily 
done in constant time. Depending on whether p crosses the boundary of C from 
its top side, right side, or from both sides simultaneously, there are three cases. 

1. lip crosses the boundary of G from the top side and p does not cross the right 
side of C, then p enters into a new cell that is on top of G. We update C to the 
new cell. We increase the index Ip{Yi) by one, but keep Ip{Xi) unchanged. 
Clearly, the above two indices are correctly updated and Ic{Xi) = Ip{Xi) 
and IciXi) — HpXi) foi' the new cell G. Again, by Lemma [TJ the function 
Edi{x,y,C) for the new cell G can be computed in constant time. We add 
the new function to Hi{R,L). 

2. lip crosses the boundary of G from the right side and p does not cross the top 
side of C, then p enters into a new cell that is on right of G. The algorithm 
in this case is similar to the above case and we omit the discussions. 

3. The remaining case is when p crosses the boundary of G through the top 
right corner of C . In this case, p enters into the northeast neighboring cell of 
G. We first add to Hi{R, L) the supporting planes of the surface patches of 
Edi{x^y) defined on the top neighboring cell and the right neighboring cell 
of C , which can be computed in constant time as the above two cases. Then, 
we update C to the new cell p is entering. We increase each of Ip{Xi) and 
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Ip{Yi) by one. Again, the two indices are correctly updated for the new cell 

C. Finally, we compute the new function Edi(a;, y, C) and add it to Hi{R, L). 

When the algorithm stops, Hi{R,L) is computed. In general, during the 
procedure of moving p on L, we spend constant time on finding each supporting 
plane of Hi{R,L). Therefore, the total running time of the entire algorithm is 
0{ai + bi). The lemma thus follows. □ 

With the preceding lemma, we have the following result. 

Theorem 1. The decision problem can be solved in 0(J27=ii^i + ^0) 

4 Computing the Rectilinear Center 

In this section, with the help of our decision algorithm in Section [Sj we compute 
the rectilinear center q* in 0{mn) time. 

As discussed in Section o our algorithm is a prune-and-search algorithm 
that has O(logn) “outer” recursive steps each of which has 0{logm) “inner” re¬ 
cursive steps. In each outer recursive step, the algorithm prunes at least \V\/32 
uncertain points of V such that these uncertain points are not relevant for com¬ 
puting q*. After O(logn) outer recursive steps, there will be only a constant 
number of uncertain points remaining in V. Each outer recursive step runs in 
0{m\V\) time, where |7^| is the number of uncertain points remaining in V. In 
this way, the total running time of the algorithm is Oimn). 

Each outer recursive step is another recursive prune-and-search algorithm, 
which consists of 2 -|- logm inner recursive steps. Let X = and y = 

Hence, \X\ = = mn. We maintain a rectangle R = [xi,X 2 ',yi,y 2 ] 

that contains q*. Initially, R is the entire plane. In each inner recursive step, we 
shrink R such that the x-range [xi,X 2 ] (resp., y-range [j/i, 2 / 2 ]) of the new R only 
contains half of the values of X (resp., 3^) in the a;-range (resp., y-range) of the 
previous R. In this way, after logm -I- 2 inner recursive steps, the a;-range (resp., 
y-range) of R only contains at most n/4 values of X (resp., y). At this moment, 
a key observation is that there is a subset V* of at least n/2 uncertain points, 
such that for each Pi € V*, R is contained in the interior of a cell of the grid Gi, 
i.e., the x-range (resp., y-range) of R does not contain any value of Xi (resp., 
Yi). Due to the observation, we can use a pruning procedure similar to that in 
m to prune at least |'P*|/16 > n/32 uncertain points. 

In the following, in Section 14.11 we give our algorithm on pruning the values 
of X and y to obtain V*. In Section 14.21 we prune uncertain points of V *. 

4.1 Pruning the Coordinate Values of X and y 

Consider a general step of the algorithm where we are about to perform the j-th 
inner recursive step for I < y < logm -I- 2. Our algorithm maintains the following 
algorithm invariants. (I) We have a rectangle R^~^ = y^~^,y^~^] 

that contains q*. (2) For each I < i < n, the index /^j-i(Ai) of the predecessor 
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of x{ ^ in Xi is known, and so is the index lyj-^ (Yi). (3) We have a sublist X/ ^ 
of Xi that consists of all values of Xi in and a sublist of Yi 

that consists of all values of Yi in [yi~^,y 2 ~^]- Note that these sublists can be 
empty. (4) < mn/2^~^ and < mn/2^~^, where 

and = Uf^^Yf-\ 

Initially, we set = [— 00 , + 00 ; — 00 , +cx)], X^ = Xi and = Yi for each 
1 < i < n, with <41° = X and = y. It is easy to see that before we start the 
first inner recursive step for j = 1, all the algorithm invariants hold. 

In the sequel, we give the details of the j-th inner recursive step. We will 
show that its running time is 0(mn/2^ +n) and all algorithm invariants are still 
maintained after the step. 

Let Xm be the median of X^~^ and ym be the median of y^~^. Both Xm and 
ym can be found in 0{\X^~^ \ + |W“^|) time. 

For each 1 < i < n, let = I{x{~^) — + 1 and = I{y 2 ~^) — 

^{y{~^) + 1 - Observe that = \Xf~^ \ + 1 and 6^”^ = + I. 

Let X* and y* be the x- and y-coordinates of q*, respectively. 

We first determine whether x* > Xm, x* < Xm, or x* = Xm- This can 
be done by applying our decision algorithm on and L with L being the 

vertical line x = Xm- By Theorem [U the running time of our decision algorithm 
is which is 0(n + \X^~^\ + 

Note that if x* = Xm, then according to our decision algorithm, q* will be 
found by the decision algorithm and we can terminate the entire algorithm. Oth¬ 
erwise, without loss of generality, we assume x* > Xm- We proceed to determine 
whether y* > y^ or y* < ym, or y* = ym by applying our decision algorithm on 
and L with L being the horizontal line y = ym- Similarly, if y* = ym, then 
the decision algorithm will find q* and we are done. Otherwise, without loss of 
generality we assume y* > ym- The above calls our decision algorithm twice, 
which takes 0{n X -I- time in total. 

Now we know that q* is in the rectangle [xm,ym, 2 / 2 ”^]- We let W = 
[xj,X 2 ;yi,y 2 ] be the above rectangle, i.e., x\ = Xm, x ^2 = ^ 2 ^ , vi = ym, and 
yi ~ yi~ ■ Clearly, the first algorithm invariant is maintained. 

We further proceed as follows to maintain the other three invariants. 

For each 1 < * < n, by scanning the sorted list X\~'^, we compute the 
index of the predecessor of x\ in Xi (each element of Xf~^ maintains its 

original index in X^), and similarly, by scanning the sorted list Y^~^, we compute 
the index I^j (Yi). Computing these indices in all Xi and Yi for i = I, 2,..., n 
can be done in 0{\X^~^ \ + time. This maintains the second algorithm 

invariant. 

Next, for each 1 < f < n, we scan to compute a sublist Xf, which 

consists of all values of Xf~ in [x{, X 2 ], and similarly, we scan Y^~ to compute a 
sublist Y-^, which consists of all values of Y^^~^ in [y( , y^]- Computing the lists X/ 
and Yf for alH = I, 2,..., n as above can be done in overall 0{\X^~^ \ + |) 

time. This maintains the third algorithm invariant. 
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Let ~ ■ According to our above algo¬ 
rithm, \X^\ < \X^~'^\/2 and 13^-^| < \y^~^\/2. Since < nm/2'^“^ and 

< nm/2-'“^, we obtain \X^\ < nm/2-' and l^-^l < nm/2K Hence, the 
fourth algorithm invariant is maintained. 

In summary, after the j-th inner recursive step, all four algorithm invariants 
are maintained. Our above analysis also shows that the total running time is 
0{n + -I- which is 0(nml2^ X n). 

We stop the algorithm after the t-th inner recursive step, for t = 2 -|- logm. 
The total time for all t steps is thus 0(X]j=i(^ + 1 ^ 71 / 2 ^)) = 0{mn). 

After the t-th step, by our algorithm invariants, the rectangle i?* contains 
q\ and |A'‘| < mn/2^ = n/4 and |y | < mn/2^ = n/4. 

We say that an uncertain point Pi is prunable if both Xf and Y* are empty 
(and thus R* is contained in the interior of a cell of Gi). Let P* denote the 
set of all prunable uncertain points of V. The following is an easy but crucial 
observation. 

Observation 2 \V*\ > n/2. 

Proof. Since X* < n/4, among the n sets Xf for i = 1, 2,..., n, at most n/4 
of them are non-empty. Similarly, since y* < n/4:, among the n sets Y/ for 
i = 1, 2,..., n, at most n/4 of them are non-empty. Therefore, there are at most 
n/2 uncertain points Pi € P such that either X/ or Y/ is non-empty. This implies 
that there are at least n/2 prunable uncertain points in P. □ 

After the t-th inner recursive step, the set P* can be obtained in 0(n) time 
by checking all sets Xj and Y/ for i = 1,2,... ,n and see whether they are empty. 

The reason we are interested in prunable uncertain points is that for each 
prunable uncertain point Pi of P*, since i?* contains q* and Rt is contained in 
a cell Ci of Gi, there is only one surface patch of Edi(x, y) (i.e., the one defined 
on Ci) that is relevant for computing q*. Let hi denote the supporting plane of 
the above surface patch. We call hi the relevant plane of Pi. Note that we can 
obtain hi in constant time. Indeed, observe that the predecessor index Ic^iXi) 
is exactly I^t{Xi), which is known by our algorithm invariants. Similarly, the 
index Ici(Yi) is also known. By Lemma [U the function Edi{x,y,Ci), which is 
also the function of hi, can be obtained in constant time. Hence, the relevant 
planes of all prunable uncertain points of P* can be obtained in 0(n) time. 


Remark. One may wonder why we did not perform the inner recursive steps 
for t = log mn times (instead of t = 2 -|- log m time) so that X* and y* would 
each have a constant number of values in the range of R. The reason is that 
based on our analysis, that would take 0{mn + nlognm) time, which may not 
be bounded by 0{mn) (e.g., when m = o(logn)). In fact, performing the inner 
recursive steps for t = 2 -t- log m times such that X* and y* each have at most ^ 
values in the range of R is an interesting and crucial ingredient of our techniques. 


12 


4.2 Pruning Uncertain Points from V* 

Consider a prunable uncertain point Pi of V*. Recall that Hi is the set of sup¬ 
porting planes of all surface patches of Edi(a;, y). The above analysis shows that 
among all planes in Hi, only the relevant plane hi is useful for determining q*. 
In other words, the point p*, as a lowest point of all planes in is 

also a lowest point of the planes in the union of and Up^^-pyp*Hi. This 

will allow us to prune at least |P*|/16 uncertain points from V*. The idea is 
similar to Megiddo’s pruning scheme for the 3D LP algorithm in m- 

For each Pi € V*, its relevant plane hi is also considered as a function in 
the xy-plane. Arrange all uncertain points of V* into |P*|/2 disjoint pairs. Let 
D{V*) denote the set of all these pairs. For each pair {Pi,Pj) G D{P*), if the 
value of the function hi at any point of R* is greater than or equal to that of hj , 
then Pj can be pruned immediately; otherwise, we project the intersection of hi 
and hj on the a:?/-plane to obtain a line Lij dividing Eh into two parts, such that 
hi > hj on one part and hi < hj on the other. 

Let £ denote the set of the dividing lines Lij for all pairs of D{V*). Let Lm 
be the line whose slope has the median value among the lines of £. We transform 
the coordinate system by rotating the x-axis to be parallel to Lm (the p-axis does 
not change). For ease of discussion, we assume no other lines of £ are parallel to 
Lm (the assumption can be easily lifted; see [H]). In the new coordinate system, 
half the lines of £ have negative slopes and the other half have positive slopes. 
We now arrange all lines of £ into disjoint pairs such that each pair has a line 
of a negative slope and a line of positive slope. Let D{jC) denote the set of all 
these line pairs. 

For each pair [Li,Lj) € D{C), we define yij as the ^-coordinate of the in¬ 
tersection of Li and Lj. We find the median ym of the values for all pairs in 
D{C). Let X* and y* respectively be the x- and y-coordinate of q* in the new 
coordinate system. We determine in 0(mn) time whether y* > ym, y* < Vm or 
y* = Vm by using our decision algorithm (here an Oimn) time decision algorithm 
is sufficient for our purpose). If y* = ym, then our decision algorithm will find 
q* and we can terminate the algorithm. Otherwise, without of loss generality, 
we assume y* < ym- 

Let D'{C) denote the set of all pairs {Li,Lj) of D{C) such that yij > ym- 
Note that \D'{C)\ > |£)(£)|/2. For each pair {Li,Lj) G D'(£), let Xij be the 
x-coordinate of the intersection of Li and Lj. We find the median Xm of all such 
Xij’s. By using our decision algorithm, we can determine in Oimn) time whether 
X* > Xm, X* < Xm, Or X* = Xm- K X* = Xm, our decision algorithm will find q* 
and we are done. Otherwise, without loss of generality, we assume x* < Xm- 

Now for each pair [Li,Lj) of D'{£) with x^ > Xm and yij > ym (there are 
at least |£>'(£)|/2 such pairs), we can prune either Pi or Pj, as follows. Indeed, 
one of the lines in such a pair [Li,Lj), say Li, has a negative slope and does not 
intersect the region R = {(x, y) | x < Xm, y < ym} (e.g., see Fig. [5]). Suppose Li 
is the dividing line of two relevant planes h^^ and hk 2 of two uncertain points 
Pfej and Pk 2 of V*. It follows that either hk^ > hk 2 or hk^ < holds on the 
region R. Since q* G R, one of Pk^ and Pk 2 can be pruned. 
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Figure 2. The intersection of Li and Lj is in the first quarter of the intersection of 
X = Xm and y = ym while q* is in the interior of the third quarter. 

As a summary, the above pruning algorithm prunes at least |7^*|/16 > n/32 
uncertain points and the total time is 0{mn). 

4.3 Wrapping Things Up 

The algorithm in the above two subsections either computes q* or prunes at least 
n/32 uncertain points from V in 0{mn) time. We assume the latter case hap¬ 
pens. Then we apply the same algorithm recursively on the remaining uncertain 
points for 0(log n) steps, after which only a constant number of uncertain points 
remain. The total running time can be described by the following recurrence: 
T(m, n) = T{m, + 0(mn). Solving the recurrence gives T{m, n) = 0{mn). 

Let V' be the set of the remaining uncertain points, with \V'\ = 0{1). Hence, 
the rectilinear center q* is determined by V'. In other words, q* is also a recti¬ 
linear center of V'. In fact, like other standard prune-and-search algorithms, the 
way we prune uncertain points of V guarantees that any rectilinear center of V 
is also a rectilinear center of V' ^ and vice versa. By using an approach similar to 
that in Section ITTI Lemma |3] finally computes q* based on V' in 0(m) time. 

Lemma 3. The rectilinear center q* can be computed in 0{m) time. 

Proof. Let c = \V'\, which is a constant. Let X' = and y' = Up.gp'U. 

We apply the same recursive algorithm in Section |4T] on X' and y' for O(logm) 
steps, after which we will obtain a rectangle R such that R contains q* and for 
each Pi G V, the a:-range (resp., y-range) of R only contains a constant number 
of values of Xi (resp., U), and thus R intersects a set Gi{R) of only a constant 
number of cells of Gi. Therefore, for each Pi G V, only the surface patches 
of Edi{x,y) defined on the cells of Gi{R) are relevant for computing q*. The 
supporting planes of these surface patches can be determined immediately after 
the above O(logm) recursive steps. By the same analysis as in Section HlTl all 
above can be done in 0{c ■ m) time. 

The above found 0(c) “relevant” supporting planes such that q* corresponds 
to a lowest point in the upper envelope of them. Consequently, q* can be found 
in 0(c) time by applying the linear-time algorithm for the 3D LP problem [15] 
on these 0(c) relevant supporting planes. □ 

This finishes our algorithm for computing g*, which runs in 0{mn) time. 
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Theorem 2. A rectilinear center q* of the uncertain points of V in the plane 
can he computed in 0{mn) time. 

5 Concluding Remarks 

In this paper, we refine the prune-and-search technique |15j to solve in linear 
time the rectilinear one-center problem on uncertain points in the plane. Note 
that the problem can also be considered as the one-center problem on uncertain 
points in the plane under the Li distance metric. Since the Loo and Li metrics 
are closely related to each other (by rotating the coordinate axes by 45°), the 
same problem under the Loo metric can be solved in linear time as well. 

The Euclidean version of the problem seems more natural. Unfortunately, 
even if V contains only one uncertain point Pi and all locations of Pi have the 
same probability, finding a center for Pi is essentially the 1-median problem in 
the plane, which is known as the Weber problem and no exact algorithm exists 
for it due to the computation challenge [3]. 
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